% This code examines the main results with different nu.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc

addpath(genpath('Parameter'))
addpath(genpath('Utility'))

addpath(genpath('jplv7'))

%% ** Import data
Data_sum_file = 'Data/Data_Sum.xlsx';
Data_sum_m = readtable(Data_sum_file,'ReadVariableNames',1,'sheet','Monthly','basic',true);
Data_sum_q = readtable(Data_sum_file,'ReadVariableNames',1,'sheet','Quarterly','basic',true);
Data_sum_a = readtable(Data_sum_file,'ReadVariableNames',1,'sheet','Annual','basic',true);
Data_sum_pred = readtable(Data_sum_file,'ReadVariableNames',1,'sheet','Pred','basic',true);

dates = gen_eomdate(Data_sum_q.yyyy(1) , Data_sum_q.yyyy(end) , 1 , 12);%1871m1 to 2019m12

% Some variables to be used
dq_g = Data_sum_q.dg;
realltg = Data_sum_q.realltg;
realmtg = Data_sum_q.realmtg;
exretexp = Data_sum_q.exretexp;
l_q = size(Data_sum_q , 1);

nu_range = (0.005 : 0.001 : 0.2);

%% ** Store result
% Table 1
res_table_1 = nan(3 , length(nu_range));%First row coefficient, second row p-value. Using i.i.d. bootstrap

%res_table_1_2 = nan(2 , length(nu_range));%First row coefficient, second row p-value. Using simple EWC without adjustment

res_table_2 = nan(2 , length(nu_range));%First row coefficient, second row p-value

res_table_3 = nan(2 , length(nu_range));%First row coefficient, second row p-value

%% ** Redo analysis
for i = 1 : length(nu_range)
    %% * Construct new experienced dividend growth
    nu_temp = nu_range(i);
    
    expd_temp = [nan(nonnan_first(dq_g) - 1 , 1);...
                 expweightedavg(dq_g(nonnan_first(dq_g) : end) , nu_temp)];
             
    %% * Table 1
    % Only use single predictor
    Nboot = 20000;
    const = 1;
    seweight = 1;%NW
    block_bs = 1;%1 for i.i.d., 2 for cicular, 3 for stationary

    l_q_1 = 372;
    nwlags_1 = round(1.3 * (l_q_1 - 1)^(1/2));

    exret_1 = Data_sum_pred.exret(end - l_q_1 : end);
    res_temp_1 = pred_VARboot(exret_1 , expd_temp(end - l_q_1 : end) , const , nwlags_1 , seweight , Nboot , block_bs);
    
    res_table_1(1 , i) = res_temp_1.badj(2);
    res_table_1(2 , i) = res_temp_1.p2null(2);
    res_table_1(3 , i) = res_temp_1.b(2);
    

    %% * Table 2
    % Use only mu_d
    LHS_2A = exretexp(2 : end);%Lagged one period to match timing
    RHS_2A = expd_temp(1 : end - 1);
    
    % Estimate
    [b_temp_2A , se_temp_2A , p_temp_2A , ~] = olsnan_EWC(LHS_2A , RHS_2A , 1);
    res_table_2(1 , i) = b_temp_2A(2);
    res_table_2(2 , i) = p_temp_2A(2);
    
    %% * Table 3
    LHS_3 = realltg(2 : end);
    RHS_3 = expd_temp(1 : end - 1);
    control_3 = realmtg(2 : end);
    
    % Estimate
    [b_temp_3 , se_temp_3 , p_temp_3 , ~] = olsnan_EWC(LHS_3 , [RHS_3 , control_3] , 1);
    res_table_3(1 , i) = b_temp_3(2);
    res_table_3(2 , i) = p_temp_3(2);
    
end

%% ** Visualization
[~ , pos_ori] = min(abs(nu_range - 0.018));

[~ , pos_sig] = min(res_table_1(2 , :));

disp('nu with highest significance')
disp(nu_range(pos_sig))

%% * Table 1
figure
yyaxis left
h_1 = plot(nu_range(1 : end) , res_table_1(1 , 1 : end) , 'LineWidth' , 2); hold on
h_2 = plot(nu_range(pos_ori) , res_table_1(1 , pos_ori) , '.' , 'MarkerSize', 25);
h_2.Annotation.LegendInformation.IconDisplayStyle = 'off';
ylabel('Coeff. (bias-adjusted)')
yyaxis right
g_1 = plot(nu_range(1 : 2 : end) , smooth(res_table_1(2 , 1 : 2 : end)) , 'LineWidth' , 2); hold on
g_2 = plot(nu_range(pos_ori) , res_table_1(2 , pos_ori) , '.' , 'MarkerSize', 25);
g_2.Annotation.LegendInformation.IconDisplayStyle = 'off';
ylabel('Bootstrapped p-value')
yline(0.05 , '--r')
xlabel('Gain parameter')
legend('Estimate' , 'p-value' , 'Location' , 'northwest')
set(gca,'TickLabelInterpreter', 'latex')
saveas(gcf,'Figures/Robust_Table1','epsc')
savetightfigure('Figures/Robust_Table1')

%% * Table 2
figure
yyaxis left
h_1 = plot(nu_range(1 : end) , res_table_2(1 , 1 : end) , 'LineWidth' , 2); hold on
h_2 = plot(nu_range(pos_ori) , res_table_2(1 , pos_ori) , '.' , 'MarkerSize', 25);
h_2.Annotation.LegendInformation.IconDisplayStyle = 'off';
ylabel('Coeff.')
yyaxis right
g_1 = plot(nu_range(1 : 2 : end) , smooth(res_table_2(2 , 1 : 2 : end)) , 'LineWidth' , 2); hold on
g_2 = plot(nu_range(pos_ori) , res_table_2(2 , pos_ori) , '.' , 'MarkerSize', 25);
g_2.Annotation.LegendInformation.IconDisplayStyle = 'off';
ylabel('EWC p-value')
yline(0.05 , '--r')
xlabel('Gain parameter')
legend('Estimate' , 'p-value' , 'Location' , 'northeast')
set(gca,'TickLabelInterpreter', 'latex')
saveas(gcf,'Figures/Robust_Table2','epsc')
savetightfigure('Figures/Robust_Table2')


%% * Table 3
figure
yyaxis left
h_1 = plot(nu_range , res_table_3(1 , :) , 'LineWidth' , 2); hold on
h_2 = plot(nu_range(pos_ori) , res_table_3(1 , pos_ori) , '.' , 'MarkerSize', 25);
h_2.Annotation.LegendInformation.IconDisplayStyle = 'off';
ylabel('Coeff.')
yyaxis right
g_1 = plot(nu_range , res_table_3(2 , :) , 'LineWidth' , 2); hold on
g_2 = plot(nu_range(pos_ori) , res_table_3(2 , pos_ori) , '.' , 'MarkerSize', 25);
g_2.Annotation.LegendInformation.IconDisplayStyle = 'off';
ylabel('EWC p-value')
yline(0.05 , '--r')
xlabel('Gain parameter')
legend('Estimate' , 'p-value' , 'Location' , 'northeast')
set(gca,'TickLabelInterpreter', 'latex')
saveas(gcf,'Figures/Robust_Table3','epsc')
savetightfigure('Figures/Robust_Table3')

